The ISS operates within the LEO environment and, given its size, at-mospheric drag results in noticeable orbit decay that requires periodic corrections via thrusters. Toexamine this in more detail, suppose that the ISS is at an altitude of 350 km and in a circular orbit.If the ballistic coefficient for the ISS is estimated as B= 165 kg/m^2, estimate the orbit decay over aperiod of six months.
#Rex Calabrese
import numpy as np
import numpy.linalg as la
import scipy as sp
from scipy import stats
from scipy.integrate import odeint
from scipy.integrate import solve_ivp
from scipy.signal import find_peaks
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from IPython.display import Image
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
from matplotlib import rc
rc('text', usetex=False)
fs_default = (11,4)
Image("images/hw8_diagram_1.png",width=600)
from IPython.display import Video
Video('Drag_animation_2_nvidiah264.mp4',embed=True,width=600)
me = 5.974e24 #kg
G = 6.67259e-20 #km^3 · kg^-1 · s^-2
re = 6378.0 #km
mu = G*me # km^3 · s^-2
A = 1. * (1. / 1000)**2 #km
Cd = 2. #unitless
m = 100 #kg
#Ballistic Coefficent
B = m / (Cd * A) #kg / km^2
print('Ballistic Coefficent: B = {0:0.3e} kg / km^2'.format(B))
r0 = [5873.40, -658.52, 3007.49] #km
v0 = [-2.90, 4.09, 6.14] #km/s
u0 = np.hstack((r0, v0))
la.norm(r0)
def rho_fn(z):
return 0
# %%
def rho_fn(z):
if z < 100:
return float(1e9)
if z >= 100:
return float(2e9 * z**(-8.284))
#Density Model: Interpolation Method
#Load data from MSISE90 model
import pandas as pd
MSISE90_pd = pd.read_excel ('density_tables/nasa_MSISE90.xlsx') #kg/km^3
MSISE90=MSISE90_pd.values
def rho_fn(z):
return np.interp(z, MSISE90[:,0], MSISE90[:,1]) / 1000 # / 1000 is a hack delete
fig00, ax00 = plt.subplots(figsize=fs_default)
for z in np.linspace(100,150,50):
ax00.scatter(z,rho_fn(z),c='grey')
plt.show()
print("SANITY CHECK")
Image("images/orbitalDrag_0.png",width=500)
print('Model says')
print('ρ(400 km) = {0:0.3e} kg / m^3'.format(rho_fn(400) * 1000**-3)) #kg/m^3
#print('... looks reasonable')
'{0:0.5e}'.format(0.0000082195)
omega_earth = [0., 0., 7.27e-5] #rad/s
def v_rel_fn(r,v,omega_earth):
return v - np.cross(omega_earth,r)
def a_drag_fn(r,r_mag,v):
v_rel = v - np.cross(omega_earth,r)
v_rel_mag = la.norm(v_rel)
rho = rho_fn(r_mag - re)
return -0.5 * (rho / B) * v_rel * v_rel_mag
# %% State Vector Derivative Function
def deriv(u, dt):
r = u[0:3]
r_mag = la.norm(r)
v = u[3:6]
n = -mu / r_mag**3
a_d = a_drag_fn(r,r_mag,v)
return [u[3], # dotu[0] = u[3]'
u[4], # dotu[1] = u[4]'
u[5], # dotu[2] = u[5]'
u[0] * n + a_d[0], # dotu[3] = u[0] * n
u[1] * n + a_d[1], # dotu[4] = u[1] * n
u[2] * n + a_d[2]] # dotu[5] = u[2] * n
def twobody(dt,u0):
u = odeint(deriv, u0, dt, atol=1.0E-8) #1E-8
r = u[:,0:3]
v = u[:,3:6]
return r,v
Y2S = 365 * 24 * 3600
D2S = 24*3600
Y2D = 365
#ORIGINAL: 5 years
nDays = 5 * Y2D
#SMALLER
#nDays = 1 * Y2D
ElapsedTime = nDays*D2S
TimePeriod = int(ElapsedTime / 5)
CALCULATE_HIGH_RESOLUTION = True
if CALCULATE_HIGH_RESOLUTION == True:
dt_hires = np.arange(0.0, ElapsedTime, 1)
r1_hires = np.zeros((dt_hires.size,3))
v1_hires = np.zeros((dt_hires.size,3))
u0_i = u0
for i in range(5): #5
i_L = i * TimePeriod
i_R = (i+1) * TimePeriod
dt_i = dt_hires[i_L:i_R]
r1_hires[i_L:i_R],v1_hires[i_L:i_R] = twobody(dt_i,u0_i)
u0_i = np.hstack((r1_hires[i_R-1,:],v1_hires[i_R-1,:]))
LOAD_HIGH_RESOLUTION = False #127 million time steps
if LOAD_HIGH_RESOLUTION == True:
r1_hires = np.load('./results_export/full/2/r1.npy')
v1_hires = np.load('./results_export/full/2/v1.npy')
dt1_hires = np.load('./results_export/full/2/dt.npy')
r1_hires[-1]
dt_hires[-1]
#fix this for solving vs loading data
dt1_hires = dt_hires
# %% Chop off the inside of the high resolution trajectory data
ends = 10000
r1_ends = np.vstack((r1_hires[0:ends,:],r1_hires[-ends:,:]))
v1_ends = np.vstack((v1_hires[0:ends,:],v1_hires[-ends:,:]))
dt1_ends = np.vstack((dt1_hires[0:ends],dt1_hires[-ends:]))
NU = np.zeros(dt1_ends.size)
# %% Keep the tail end of hires trajectory data.
r1 = r1_hires[-100000:] #100k
v1 = v1_hires[-100000:] #100k
dt1 = dt1_hires[-100000:] #100k
dt1_days = dt1 / D2S
r1_mag = la.norm(r1,axis=1)
v1_mag = la.norm(v1,axis=1)
# %% Downsample hires trajectory data.
sample_interval = int(D2S / 10)
r2 = r1_hires[0::sample_interval,:]
del r1_hires
v2 = v1_hires[0::sample_interval,:]
del v1_hires
dt2_days = dt1_hires[0::sample_interval] / D2S
del dt1_hires
dt2 = dt2_days / Y2D
# Subscript 2 data is downsampled.
# Subscript 1 data is clipped.
h2 = np.cross(r2,v2,axis=1)
r2_mag = np.linalg.norm(r2, axis=1)
v2_mag = np.linalg.norm(v2, axis=1)
eccvec2 = np.cross(v2, h2)/mu - np.divide(r2,r2_mag[:,None])
eccmag2 = la.norm(eccvec2,axis=1)
energy2 = (v2_mag**2 / 2) - (mu/r2_mag)
a2 = -mu/(2*energy2)
r2_a = a2 * (1 + eccmag2)
r2_p = a2 * (1 - eccmag2)
fig0, ax0 = plt.subplots(figsize=fs_default)
ax0.plot(dt2,a2,c='grey')
slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,a2)
line = slope*dt2+intercept
ax0.plot(dt2,line,'r--',lw=1)
ax0.set_title('Semi-Major axis : $\Delta a$/$\Delta t$={0:0.2f} km/yr'.format(slope,intercept))
ax0.set_xlabel('Time (yr)')
ax0.set_ylabel('Semi-Major Axis (km)')
ax0.legend(['e','linear regression'])
ax0.grid(True)
plt.show()
dt2[-1]
print('semi major axis decay over {0:0.2f} years'.format(dt2[-1] ))
print('Δa = {0:0.3f} km'.format(a2[-1] - a2[0] ))
fig1, ax1 = plt.subplots(figsize=fs_default)
ax1.plot(dt2,eccmag2,c='grey')
slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,eccmag2)
line = slope*dt2+intercept
ax1.plot(dt2,line,'r--',lw=1)
ax1.set_title('Eccentricity : $\Delta e$/$\Delta t$={0:0.2e} yr$^-1$'.format(slope,intercept))
ax1.set_xlabel('Time (yr)')
ax1.set_ylabel('Eccentricity (km/km)')
ax1.legend(['e','linear regression'])
plt.show()
fig3, ax3 = plt.subplots(figsize=fs_default)
c_rp = 'orange'
c_ra = 'purple'
#apoapsis
ax3.plot(dt2,r2_a - re,c=c_ra,lw=3,label='$apogee$')
slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,r2_a - re)
line = slope*dt2+intercept
ax3.plot(dt2,line,'r--',lw=1,label='$\Delta r_a$/$\Delta t$ = {0:0.1f} km/yr'.format(slope))
#perigee
ax3.plot(dt2,r2_p - re,c=c_rp,lw=3,label='$perigee$')
slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,r2_p - re)
line = slope*dt2+intercept
ax3.plot(dt2,line,'b--',lw=1,label='$\Delta r_p$/$\Delta t$ = {0:0.1f} km/yr'.format(slope))
ax3.set_title('Atmospheric Drag: Perigee and Apoapsis Values')
ax3.set_ylabel('$Altitude$ (km)')
ax3.set_xlabel('Time (years)')
ax3.legend()
ax3.grid(True)
plt.show()
# %% Plot flat orbit with r_a and r_p
r2_a_shift = r2_a-r2_a[0]
r2_p_shift = r2_p-r2_p[0]
fig4, ax4 = plt.subplots(figsize=fs_default)
#apoapsis
ax4.plot(dt2,r2_a_shift,c=c_ra,lw=3,label='$r_a$')
slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,r2_a_shift)
line = slope*dt2+intercept
ax4.plot(dt2,line,'r--',lw=1,label='$\Delta r_a$/$\Delta t$ = {0:0.1f} km/yr'.format(slope))
#perigee
ax4.plot(dt2,r2_p_shift,c=c_rp,lw=3,label='$r_p$')
slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,r2_p_shift)
line = slope*dt2+intercept
ax4.plot(dt2,line,'b--',lw=1,label='$\Delta r_p$/$\Delta t$ = {0:0.1f} km/yr'.format(slope))
ax4.set_title('Atmospheric Drag: Perigee and Apoapsis Values (arbitrary initial values)')
ax4.set_ylabel('$r$ (km)')
ax4.set_xlabel('Time (years)')
ax4.legend()
plt.show()
The above plot shows a relatively steady perigee. We know the drag force is largest at perigee due to increased atmospheric density, similar to an instantaneous ΔV maneuver at perigee, this location remains 'pinned'.
# %% 'Micro View': FIND r_a and r_p locations, and plot over r vs t.
r1_apo = find_peaks(r1_mag)
r1_peri = find_peaks(-r1_mag)
plt.figure(figsize=fs_default)
#Plot r
plt.plot(dt1_days,r1_mag,color='k')
#Plot the found r_a,t_p locations
for i in range(r1_peri[0].size - 1): #
j = r1_peri[0][i]
j1 = r1_apo[0][i]
plt.scatter(dt1_days[j],r1_mag[j], c=c_ra)
plt.scatter(dt1_days[j1],r1_mag[j1], c=c_rp)
plt.legend(['|r|','Apogee','Perigee'],loc=1)
plt.title('Test: Find Perigee and Apoasis Values Through Peak Searching')
plt.show()
# %% RECALCULATE ACCEL (#Using r1, v1 ,dt1 )
a1 = np.zeros((dt1.size,3))
for i in range(dt1.size):
a1[i,:] = a_drag_fn(r1[i,:],r1_mag[i],v1[i,:])
a1_mag = la.norm(a1,axis=1)
a1[1]
# Recalculate V_rel
v1rel = np.zeros((dt1.size,3))
for i in range(dt1.size):
v1rel[i,:] = v_rel_fn(r1[i,:],v1[i,:],omega_earth)
v1rel_mag = la.norm(v1rel,axis=1)
dt1_yr = dt1 / Y2S
dt1_day = dt1 / D2S
fig5, ax5 = plt.subplots(3, figsize=(11,8), sharex=True)
ax5[0].plot(dt1_day, r1[:,0],c='r',label='$r_x$')
ax5[0].plot(dt1_day, r1[:,1],c='g', label='$r_y$')
ax5[0].plot(dt1_day, r1[:,2],c='b', label='$r_z$')
ax5[0].plot(dt1_day, r1_mag[:],c='k', label='|r|')
ax5[1].plot(dt1_day, v1rel[:,0],c='r', label='$v_x$')
ax5[1].plot(dt1_day, v1rel[:,1],c='g', label='$v_y$')
ax5[1].plot(dt1_day, v1rel[:,2],c='b', label='$v_z$')
ax5[1].plot(dt1_day, v1rel_mag,c='k', label='|v|')
ax5[2].plot(dt1_day, a1[:,0],c='r', label='$a_x$')
ax5[2].plot(dt1_day, a1[:,1],c='g', label='$a_y$')
ax5[2].plot(dt1_day, a1[:,2],c='b', label='$a_z$')
ax5[2].plot(dt1_day, a1_mag,c='k', label='|a|')
#r_a and r_p
for i in range(r1_apo[0].size): #
j1 = r1_apo[0][i]
if i ==0:
ax5[0].scatter(dt1_day[j1],r1_mag[j1], c=c_rp, label='A')
ax5[1].scatter(dt1_day[j1],v1rel_mag[j1], c=c_rp,label='A')
ax5[2].scatter(dt1_day[j1],a1_mag[j1], c=c_rp,label='A')
else:
ax5[0].scatter(dt1_day[j1],r1_mag[j1], c=c_rp)
ax5[1].scatter(dt1_day[j1],v1rel_mag[j1], c=c_rp)
ax5[2].scatter(dt1_day[j1],a1_mag[j1], c=c_rp)
for i in range(r1_peri[0].size): #
j = r1_peri[0][i]
if i ==0:
ax5[0].scatter(dt1_day[j],r1_mag[j], c=c_ra, label='P')
ax5[1].scatter(dt1_day[j],v1rel_mag[j], c=c_ra, label='P')
ax5[2].scatter(dt1_day[j],a1_mag[j], c=c_ra, label='P')
else:
ax5[0].scatter(dt1_day[j],r1_mag[j], c=c_ra)
ax5[1].scatter(dt1_day[j],v1rel_mag[j], c=c_ra)
ax5[2].scatter(dt1_day[j],a1_mag[j], c=c_ra)
for i in range(3): ax5[i].legend(loc=1)
ax5[0].set_title('Position, Relative Velocity, Deceleration (final orbit)')
for i in range(3): ax5[i].set_xlim((5*Y2D - 0.18,5*Y2D))
ax5[2].set_xlabel('Time (days)')
ax5[0].set_ylabel('Position (km)')
ax5[1].set_ylabel('Relative Velocity (km/s)')
ax5[2].set_ylabel('Deceleration (km/s$^2$)')
plt.show()
The above plot shows the deceleration values are inversely proportional to altitude.
%matplotlib notebook
#fig = plt.figure(dpi=100) #
fig = plt.figure(figsize=(6,6),dpi=135)
#fig = plt.figure(dpi=100)
#ax = fig.add_subplot(111, projection = '3d')
ax = Axes3D(fig)
ax.set_xlim([-8000,8000])
ax.set_ylim([-8000,8000])
ax.set_zlim([-8000,8000])
###########################################
#Add Earth (Bottom Half)
ue1, ve1 = np.mgrid[0:2*np.pi:60j, 0:np.pi:60j]
xe=re*np.cos(ue1)*np.sin(ve1)
ye=re*np.sin(ue1)*np.sin(ve1)
ze=re*np.cos(ve1)
clip0 = int(xe.size / 10)
ax.plot_surface(xe, ye, ze, lw=0.1,alpha=0.2,color='deepskyblue',zorder = -10)
ax.plot_wireframe(xe, ye, ze, lw=0.1,alpha=0.99,color='deepskyblue',zorder = -10)
#Show Earth Axis of rotation.
ax.plot([0,0],[0,0],[-re*1.2,re*1.2],color='k',lw=4,ls='-',alpha=0.2,label='Axis of Rotation: Earth')
ax.scatter(0,0,0,color='k',alpha=0.8)
###########################################
#Plot First Full Orbit
clip1 = 0
clip2 = 7500
ax.plot(r1[clip1:clip2,0], r1[clip1:clip2,1], r1[clip1:clip2,2],'k',lw=1,zorder = 2.5,label='Initial Orbit')
###########################################
#apoapsis
c_v = 'mediumspringgreen'
c_a = 'hotpink'
i=0
j1 = r1_apo[0][i]
ax.scatter(r1[j1,0],r1[j1,1],r1[j1,2], c=c_ra,s=30) #show location of apoapsis
ax.quiver(0,0,0,r1[j1,0],r1[j1,1],r1[j1,2], color=c_ra,label='$r_{Apoapsis}$') # show vector to apoapsis
#Show velocity vector at apoapsis
scale1 = 1000
ax.quiver(r1[j1,0],r1[j1,1],r1[j1,2],v1rel[j1,0]*scale1,v1rel[j1,1]*scale1,v1rel[j1,2]*scale1, color=c_v)
#Show decel vector at apoapsis
scale2 = 1e25
ax.quiver(r1[j1,0],r1[j1,1],r1[j1,2],a1[j1,0]*scale2,a1[j1,1]*scale2,a1[j1,2]*scale2, color=c_a)
###########################################
#perigee
i=0
j1 = r1_peri[0][i]
ax.scatter(r1[j1,0],r1[j1,1],r1[j1,2], c=c_rp,s=30)#show location of perigee
ax.quiver(0,0,0,r1[j1,0],r1[j1,1],r1[j1,2], color=c_rp,label='$r_{Periapsis}$') # show vector to perigee
#Show velocity vector at perigee
ax.quiver(r1[j1,0],r1[j1,1],r1[j1,2],v1rel[j1,0]*scale1,v1rel[j1,1]*scale1,v1rel[j1,2]*scale1, color=c_v,label='$\\vec v$')
#Show deceleration vector at perigee (Scales of deceleration are not equal)
scale3 = 1e20
ax.quiver(r1[j1,0],r1[j1,1],r1[j1,2],a1[j1,0]*scale3,a1[j1,1]*scale3,a1[j1,2]*scale3, color=c_a,label='$ \\vec a$ (not to scale)')
###########################################
#Label
ax.set_xlabel('x (km)')
ax.set_ylabel('y (km)')
ax.set_zlabel('z (km)')
#Set ticks
ax.set_xticks([-5000,0,5000])
ax.set_yticks([-5000,0,5000])
ax.set_zticks([-5000,0,5000])
plt.legend()
ax.set_axis_off()
plt.show()
animate_rotate = True
if animate_rotate==True:
for angle in range(0, 360):
ax.view_init(0, angle)
fig.canvas.draw()
plt.pause(.001)
ax.view_init(-22,70)
The above figure is another representation of the inverse proportionality between deceleration and altitude.
%matplotlib notebook
#fig = plt.figure(dpi=100) #
#fig = plt.figure(dpi=100)
#Very good
fig = plt.figure(figsize=(6,6),dpi=135)
#ax = fig.add_subplot(111, projection = '3d')
ax = Axes3D(fig)
axlim = 7000
ax.set_xlim([-axlim,axlim])
ax.set_ylim([-axlim,axlim])
ax.set_zlim([-axlim,axlim])
ax.view_init(-22,70)
###########################################
#Add Earth (Bottom Half)
ue1, ve1 = np.mgrid[0:2*np.pi:60j, 0:np.pi:60j]
xe=re*np.cos(ue1)*np.sin(ve1)
ye=re*np.sin(ue1)*np.sin(ve1)
ze=re*np.cos(ve1)
clip0 = int(xe.size / 10)
ax.plot_surface(xe, ye, ze, lw=0.1,alpha=0.2,color='deepskyblue',zorder = -10)
ax.plot_wireframe(xe, ye, ze, lw=0.1,alpha=0.99,color='deepskyblue',zorder = -10)
#Show Earth Axis of rotation.
ax.plot([0,0],[0,0],[-re*1.2,re*1.2],color='k',lw=4,ls='-',alpha=0.2,label='Axis of Rotation')
ax.scatter(0,0,0,color='k',alpha=0.8)
###########################################
#Plot First Full Orbit
#clip1 = 0
#clip2 = 7500
#
###########################################
#apoapsis
c_v = 'mediumspringgreen'
c_a = 'hotpink'
#Label
ax.set_xlabel('x (km)')
ax.set_ylabel('y (km)')
ax.set_zlabel('z (km)')
#Set ticks
ax.set_xticks([-5000,0,5000])
ax.set_yticks([-5000,0,5000])
ax.set_zticks([-5000,0,5000])
#ax.set_axis_off()
def clear_lines():
quiv_accel.remove()
quiv_vel.remove()
quiv_radius.remove()
plt.show()
LW_vec = 3
#For parameter output in plot
line = -20000
ypos_0 = 0.7 #-0.1
xpos_0 =0.02
animation_range = 1000
for j in range(animation_range):
i = j * 30 # * 80 for display
#ORBIT
ax.plot((r1[i,0],r1[i+1,0]), (r1[i,1],r1[i+1,1]), (r1[i,2],r1[i+1,2]),'k',lw=1,color='k',zorder = 2.5,label='Orbit')
#RADIUS VECTOR
quiv_radius = ax.quiver(0,0,0,r1[i,0],r1[i,1],r1[i,2], color=c_rp,linewidths=LW_vec, label='Vector to sat')
#VELOCITY VECTOR
s1 = 1000
quiv_vel = ax.quiver(r1[i,0],r1[i,1],r1[i,2],v1rel[i,0]*s1,v1rel[i,1]*s1,v1rel[i,2]*s1,color=c_v,linewidths=LW_vec,label='Velocity')
#DECEL VECTOR
s2 = 1e20
quiv_accel = ax.quiver(r1[i,0],r1[i,1],r1[i,2],a1[i,0]*s2,a1[i,1]*s2,a1[i,2]*s2, color=c_a, linewidths=LW_vec,label = 'Drag')
text_below = '$altitude =$'+'{0:.0f} km \ndrag = {1:.1e} $km/s^2$'.format(r1_mag[i] - re,a1_mag[i] )
text_parameter_1 = ax.text(xpos_0, ypos_0 + 0*line ,0,text_below,
horizontalalignment='left',
verticalalignment='bottom',
size='large',
bbox=dict(facecolor='gray', alpha=0.3),
transform = ax.transAxes)
#SHOW
fig.canvas.draw()
plt.pause(0.0001)
if j == 0:
plt.legend(loc='lower right', borderaxespad=0.) #bbox_to_anchor=(1.05, 1),
#Save frames for video
frame_name = 'animation/frame_{}'.format(j)
plt.savefig(frame_name)
#REMOVE PREVIOUS
if j < animation_range-1:
quiv_accel.remove()
quiv_vel.remove()
quiv_radius.remove()
text_parameter_1.remove()
del quiv_accel
del quiv_vel
del quiv_radius
del text_parameter_1
#quiv_vel.set_alpha(0)
#quiv_radius.set_alpha(0)
#quiv_accel.set_alpha(0)
#
###########################################
plt.show()